Examples

Examples

As a quick hands-on introduction, we demonstrate the usage of the CoherentStructures package on some classic flow problems. For references to the original works in which the methods were developed see the respective help page.

Rotating Double Gyre

FEM-Based Methods

The following code-snippet shows how these methods can be used.

using CoherentStructures,Arpack
LL = [0.0,0.0]; UR = [1.0,1.0];
ctx = regularTriangularGrid((50,50),LL,UR)

A = x-> mean_diff_tensor(rot_double_gyre,x,[0.0,1.0], 1.e-10,tolerance= 1.e-4)
K = assembleStiffnessMatrix(ctx,A)
M = assembleMassMatrix(ctx)
λ, v = eigs(-K,M,which=:SM);

This velocity field is given by the rot_double_gyre function. The second argument to mean_diff_tensor are the times at which we average the pullback diffusion tensors. The third parameter is the step size δ used for the finite-difference scheme, tolerance is passed to the ODE solver from DifferentialEquations.jl. In the above, A(x) approximates the mean diffusion tensor given by

\[A(x) = \sum_{t \in \mathcal T}(D\Phi^t(x))^{-1} (D\Phi^t x)^{-T}\]

The eigenfunctions saved in v approximate those of $\Delta^{dyn}$

import Plots
res = [plot_u(ctx, v[:,i],colorbar=:none,clim=(-3,3)) for i in 1:6];
Plots.plot(res...,margin=-10Plots.px)
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

Looking at the spectrum, there appears a gap after the third eigenvalue:

Plots.scatter(1:6, real.(λ))
1 2 3 4 5 6 0 100 200 300 400 500 y1

We can use the Clustering.jl package to compute coherent structures from the first two nontrivial eigenfunctions:

using Clustering
numclusters=2
res = kmeans(permutedims(v[:,2:numclusters+1]),numclusters+1)
u = kmeansresult2LCS(res)
Plots.plot([plot_u(ctx,u[:,i],200,200,color=:viridis) for i in [1,2,3]]...)
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Geodesic vortices

The rotating double gyre model was introduced by Mosovsky & Meiss. It can be derived from the stream function $ \psi(x,y,t)=(1−s(t))\psiP +s(t)\psiF$ $\psi_P (x, y) = \sin(2\pi x) \sin(\pi y)$ $\psi_F (x, y) = \sin(\pi x) \sin(2\pi y)$

where $s$ is (usually taken to be) a cubic interpolating function satisfying $s(0) = 0$ and $s(1) = 1$. It therefore interpolates two double gyre flow fields, from horizontally to vertically arranged counter-rotating gyres. The corresponding velocity field is provided by the package and callable as rot_double_gyre.

Here, we demonstrate how to calculate black-hole vortices, see Geodesic elliptic material vortices for references and details.

using CoherentStructures
import Tensors, OrdinaryDiffEq, Plots

const q = 51
const tspan = collect(range(0.,stop=1.,length=q))
ny = 101
nx = 101
N = nx*ny
xmin, xmax, ymin, ymax = 0.0, 1.0, 0.0, 1.0
xspan, yspan = range(xmin,stop=xmax,length=nx), range(ymin,stop=ymax,length=ny)
P = vcat.(xspan',yspan)
const δ = 1.e-6
mCG_tensor = u -> CG_tensor(rot_double_gyre,u,tspan,δ,tolerance=1e-6,solver=OrdinaryDiffEq.Tsit5())

C = map(mCG_tensor,P)

LCSparams = (.1, 0.5, 0.01, 0.2, 0.3, 60)
vals, signs, orbits = ellipticLCS(C,xspan,yspan,LCSparams);
Detected 17 singularity candidates...
Determined 9 nondegenerate singularities...
Defined 2 Poincaré sections...
Found 2 vortices.

The results are then visualized as follows.

using Statistics
λ₁, λ₂, ξ₁, ξ₂, traceT, detT = tensor_invariants(C)
# damp "outliers"
l₁ = min.(λ₁,quantile(λ₁[:],0.999))
l₁ = max.(λ₁,1e-2)
l₂ = min.(λ₂,quantile(λ₂[:],0.995))

fig = Plots.heatmap(xspan,yspan,log10.(l₂),aspect_ratio=1,color=:viridis,
            title="FTLE-field and transport barriers",xlims=(xmin, xmax),ylims=(ymin, ymax),leg=true)
for i in eachindex(orbits)
    Plots.plot!(orbits[i][1,:],orbits[i][2,:],w=3,label="T = $(round(vals[i],digits=2))")
end
Plots.plot(fig)
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 FTLE-field and transport barriers 2.5 5.0 7.5 10.0 T = 1.0 T = 0.8

Bickley Jet

The Bickley jet flow is a kinematic idealized model of a meandering zonal jet flanked above and below by counterrotating vortices. It was introduced by Rypina et al.; cf. also del‐Castillo‐Negrete and Morrison.

The Bickley jet is described by a time-dependent velocity field arising from a stream-function. The corresponding velocity field is provided by the package and callable as bickleyJet.

Geodesic vortices

Here we briefly demonstrate how to find material barriers to diffusive transport; see Geodesic elliptic material vortices for references and details.

using CoherentStructures, Arpack
import Tensors, OrdinaryDiffEq

const q = 81
const tspan = collect(range(0.,stop=3456000.,length=q))
ny = 120
nx = div(ny*24,6)
N = nx*ny
xmin, xmax, ymin, ymax = 0.0 - 2.0, 6.371π + 2.0, -3.0, 3.0
xspan, yspan = range(xmin,stop=xmax,length=nx), range(ymin,stop=ymax,length=ny)
P = vcat.(xspan',yspan)
const δ = 1.e-6
const DiffTensor = Tensors.SymmetricTensor{2,2}([2., 0., 1/2])
mCG_tensor = u -> av_weighted_CG_tensor(bickleyJet,u,tspan,δ,
          D=DiffTensor,tolerance=1e-6,solver=OrdinaryDiffEq.Tsit5())

C̅ = map(mCG_tensor,P)

LCSparams = (.1, 0.5, 0.04, 0.5, 1.8, 60)
vals, signs, orbits = ellipticLCS(C̅,xspan,yspan,LCSparams);
Detected 318 singularity candidates...
Determined 76 nondegenerate singularities...
Defined 7 Poincaré sections...
Found 7 vortices.

The result is visualized as follows:

import Plots
using Statistics
λ₁, λ₂, ξ₁, ξ₂, traceT, detT = tensor_invariants(C̅)
l₁ = min.(λ₁,quantile(λ₁[:],0.999))
l₁ = max.(λ₁,1e-2)
l₂ = min.(λ₂,quantile(λ₂[:],0.995))
begin
    fig = Plots.heatmap(xspan,yspan,log10.(l₁.+l₂),aspect_ratio=1,color=:viridis,
            title="DBS-field and transport barriers",xlims=(0., 6.371π),ylims=(-3., 3.),leg=true)
    for i in eachindex(orbits)
        Plots.plot!(orbits[i][1,:],orbits[i][2,:],w=3,label="T = $(round(vals[i],digits=2))")
    end
end
Plots.plot(fig)
0 5 10 15 20 -3 -2 -1 0 1 2 3 DBS-field and transport barriers 2.5 5.0 7.5 10.0 T = 1.27 T = 1.27 T = 1.23 T = 1.27 T = 1.27 T = 1.28 T = 1.24

FEM-based Methods

The Bickley Jet is described by a time-dependent velocity field arising from a stream-function. Instead of using the bickleyJet function to get this velocity field, we can use the @velo_from_stream macro:

using CoherentStructures
# after this, 'bickley' will reference a Dictionary of functions
# access it via the desired signature. e.g. F = bickley[:(dU, U, p, t)]
# for the right side of the equation of variation.
bickley = @velo_from_stream stream begin
    stream = psi₀ + psi₁
    psi₀   = - U₀ * L₀ * tanh(y / L₀)
    psi₁   =   U₀ * L₀ * sech(y / L₀)^2 * re_sum_term

    re_sum_term =  Σ₁ + Σ₂ + Σ₃

    Σ₁  =  ε₁ * cos(k₁*(x - c₁*t))
    Σ₂  =  ε₂ * cos(k₂*(x - c₂*t))
    Σ₃  =  ε₃ * cos(k₃*(x - c₃*t))

    k₁ = 2/r₀      ; k₂ = 4/r₀    ; k₃ = 6/r₀

    ε₁ = 0.0075    ; ε₂ = 0.15    ; ε₃ = 0.3
    c₂ = 0.205U₀   ; c₃ = 0.461U₀ ; c₁ = c₃ + (√5-1)*(c₂-c₃)

    U₀ = 62.66e-6  ; L₀ = 1770e-3 ; r₀ = 6371e-3
end;

As we are using a periodic domain in one direction:

LL = [0.0,-3.0]; UR=[6.371π,3.0]
ctx = regularTriangularGrid((100,30),LL,UR,quadrature_order=1)
predicate = (p1,p2) -> abs(p1[2] - p2[2]) < 1e-10 && peuclidean(p1[1], p2[1], 6.371π) < 1e-10
bdata = CoherentStructures.boundaryData(ctx,predicate,[]);

Using a FEM-based method to compute coherent structures:

using Arpack,Statistics
cgfun = (x -> mean(pullback_diffusion_tensor(bickley, x,range(0.0,stop=40*3600*24,length=81),
     1.e-8,tolerance=1.e-5)))

K = assembleStiffnessMatrix(ctx,cgfun,bdata=bdata)
M = assembleMassMatrix(ctx,bdata=bdata)
λ, v = eigs(K,M,which=:SM, nev= 10)
plot_real_spectrum(λ)
2 4 6 8 10 -80 -60 -40 -20 0 y1

K-means clustering gives something we can plot:

using Clustering,Plots
n_partition = 8
res = kmeans(permutedims(v[:,2:n_partition]),n_partition)
u = kmeansresult2LCS(res)
u_combined = sum([u[:,i]*i for i in 1:n_partition])
plot_u(ctx, u_combined,200,200,bdata=bdata,
    color=:rainbow,colorbar=:none,title="$n_partition-partition of Bickley Jet")
0 5 10 15 20 -3 -2 -1 0 1 2 3 8-partition of Bickley Jet

We ran kmeans only once. To get better results, kmeans should be run several times and only the run with the lowest objective function be kept.

Geostrophic Ocean Flow

For a more realistic application, we consider an unsteady ocean surface velocity data set obtained from satellite altimetry measurements produced by SSALTO/DUACS and distributed by AVISO. The particular space-time window has been used several times in the literature.

Geodesic vortices

Here, we demonstrate how to detect material barriers to diffusive transport.

using CoherentStructures
import JLD2, OrdinaryDiffEq, Plots
###################### load and interpolate velocity data sets #############
JLD2.@load("../../examples/Ocean_geostrophic_velocity.jld2")

UI, VI = interpolateVF(Lon,Lat,Time,UT,VT)
p = (UI,VI)

############################ integration set up ################################
q = 91
t_initial = minimum(Time)
t_final = t_initial + 90
const tspan = range(t_initial,stop=t_final,length=q)
nx = 500
ny = Int(floor(0.6*nx))
N = nx*ny
xmin, xmax, ymin, ymax = -4.0, 6.0, -34.0, -28.0
xspan, yspan = range(xmin,stop=xmax,length=nx), range(ymin,stop=ymax,length=ny)
P = vcat.(xspan',yspan)
const δ = 1.e-5
mCG_tensor = u -> av_weighted_CG_tensor(interp_rhs,u,tspan,δ,p = p,tolerance=1e-6,solver=OrdinaryDiffEq.Tsit5())

C̅ = map(mCG_tensor,P)
LCSparams = (.09, 0.5, 0.05, 0.5, 1.0, 60)
vals, signs, orbits = ellipticLCS(C̅,xspan,yspan,LCSparams);
Detected 329 singularity candidates...
Determined 110 nondegenerate singularities...
Defined 8 Poincaré sections...
Found 4 vortices.

The result is visualized as follows:

using Statistics
λ₁, λ₂, ξ₁, ξ₂, traceT, detT = tensor_invariants(C̅)
l₁ = min.(λ₁,quantile(λ₁[:],0.999))
l₁ = max.(λ₁,1e-2)
l₂ = min.(λ₂,quantile(λ₂[:],0.995))
fig = Plots.heatmap(xspan,yspan,log10.(l₁.+l₂),aspect_ratio=1,color=:viridis,
            title="DBS-field and transport barriers",xlims=(xmin, xmax),ylims=(ymin, ymax),leg=true)
for i in eachindex(orbits)
    Plots.plot!(orbits[i][1,:],orbits[i][2,:],w=3,label="T = $(round(vals[i],digits=2))")
end
Plots.plot(fig)
-4 -2 0 2 4 6 -34 -33 -32 -31 -30 -29 -28 DBS-field and transport barriers 1 2 3 4 5 6 7 8 T = 1.24 T = 1.0 T = 1.19 T = 1.15

Standard Map

The "standard map" with parameter $a$ is defined on a 2-dimensional doubly 2π-periodic domain by $f(x,y) = (x+ y+ a \sin(x),y + a\cos(x))$.

For $a = 0.971635$, the standard map is implemented in CoherentStructures.standardMap, its Jacobi-matrix in CoherentStructures.DstandardMap.

See also Froyland & Junge (2015), who calculate Coherent Structures for this map.

Below are some orbits of the standard map

using Random
to_plot = []
for i in 1:50
    Random.seed!(i)
    x = rand(2)*2π
    for i in 1:500
        x = CoherentStructures.standardMap(x)
        push!(to_plot,x)
    end
end
Plots.scatter([x[1] for x in to_plot],[x[2] for x in to_plot],
    m=:pixel,ms=1,aspect_ratio=1,legend=:none)
0 1 2 3 4 5 6 0 1 2 3 4 5 6

Approximating the Dynamical Laplacian by FEM methods is straightforward:

using Tensors, Plots, Arpack, Printf
ctx = regularTriangularGrid((100,100), [0.0,0.0],[2π,2π])
pred  = (x,y) -> peuclidean(x[1],y[1],2π) < 1e-9 && peuclidean(x[2],y[2],2π) < 1e-9
bdata = boundaryData(ctx,pred) #Periodic boundary

id2 = one(Tensors.Tensor{2,2}) # 2D identity tensor
cgfun = x -> 0.5*(id2 +  Tensors.dott(Tensors.inv(CoherentStructures.DstandardMap(x))))

K = assembleStiffnessMatrix(ctx,cgfun,bdata=bdata)
M = assembleMassMatrix(ctx,lumped=false,bdata=bdata)
@time λ, v = eigs(-1*K,M,which=:SM)
Plots.plot(
    [plot_u(ctx,v[:,i],bdata=bdata,title=@sprintf("\\lambda = %.3f",λ[i]),
        clim=(-0.25,0.25),colorbar=:none)
         for i in 1:6]...)
0 1 2 3 4 5 6 0 1 2 3 4 5 6 λ = -0.000 0 1 2 3 4 5 6 0 1 2 3 4 5 6 λ = 1.103 0 1 2 3 4 5 6 0 1 2 3 4 5 6 λ = 1.103 0 1 2 3 4 5 6 0 1 2 3 4 5 6 λ = 1.501 0 1 2 3 4 5 6 0 1 2 3 4 5 6 λ = 1.501 0 1 2 3 4 5 6 0 1 2 3 4 5 6 λ = 1.689